#include <alloca.h>
#include <mml/polygon.h>
#include <mml/primitives.h>

int mml::reducePolygon(int n,
                       double* vertexData,
                       double const planeNormal[3],
                       double const& planeDistance)
{
  // Compute the distance of each vertex to the plane:
  double* distances = (double*)alloca(n * sizeof(double));
  double dmax = -DBL_MAX;
  double dmin = DBL_MAX;
  for (int i = 0; i < n; ++i) {
    distances[i] = vector_ip(3, &vertexData[3 * i], planeNormal) - planeDistance;
    if (distances[i] > dmax) {
      dmax = distances[i];
    }
    if (distances[i] < dmin) {
      dmin = distances[i];
    }
  }

  if (dmax < 0.0) {
    // Entire panel is below plane, return empty panel:
    return 0;
  }
  else if (dmin > 0.0) {
    // Entire panel is above plane and remains unchanged:
    return -1;  // -1 signifies 'no change'
  }

  // The plane intersects the panel, reduce it:
  int m = 0;
  double* newVertexData = (double*)alloca(2 * 3 * n * sizeof(double));

  // Sutherland-Hodgman algorithm.
  int currentIndex = n - 1;
  for (int nextIndex = 0; nextIndex < n; ++nextIndex) {
    if (distances[currentIndex] * distances[nextIndex] < 0.0) {
      // The current and next vertices lie on different sides of the plane.
      // Output clipped vertex (by linear interpolation):
      double factor = 1.0 / (distances[nextIndex] - distances[currentIndex]);
      int m3 = 3 * m;
      vector_reset(3, &newVertexData[m3]);
      vector_add(3, &newVertexData[m3], &vertexData[3 * currentIndex], distances[nextIndex] * factor);
      vector_add(3, &newVertexData[m3], &vertexData[3 * nextIndex], -distances[currentIndex] * factor);
      ++m;
    }

    if (distances[nextIndex] >= 0.0) {
      // Next vertex is in valid region, output it:
      vector_copy(3, &newVertexData[3 * m], &vertexData[3 * nextIndex]);
      ++m;
    }

    currentIndex = nextIndex;
  }

  vector_copy(3 * m, vertexData, newVertexData);

  return m;
}
